Majorana Fermions and Exotic Surface Andreev Bound States in Topological 
Superconductors: Application to Cu x Bi2Se3 



Timothy H. Hsieh 1 and Liang Fu 1 ' 2 

department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 
2 Department of Physics, Harvard University, Cambridge, MA 02138 

The recently discovered superconductor Cu^I^Ses is a candidate for three-dimensional time- 
reversal-invariant topological superconductors, which are predicted to have robust surface Andreev 
bound states hosting massless Majorana fermions. In this work, we analytically and numerically 
find the linearly dispersing Majorana fermions at k — 0, which smoothly evolve into a new branch 
of gapless surface Andreev bound states near the Fermi momentum. The latter is a new type 
of Andreev bound states resulting from both the nontrivial band structure and the odd-parity 
pairing symmetry. The tunneling spectra of these surface Andreev bound states agree well with a 
recent point-contact spectroscopy experiment jf! and yield additional predictions for low temperature 
tunneling and photoemission experiments. 
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The discovery of topological insulators has gener- 
ated much interest in not only understanding their 
properties and potential applications to spintronics and 
thermoelectrics but also searching for new topological 
phases. A particularly exciting avenue is topological 
superconductors [2T-fT0] . in which unconventional pairing 
symmetries lead to topologically ordered superconduct- 
ing ground states [Tllfl3] . The hallmark of a topolog- 
ical superconductor is the existence of gapless surface 
Andreev bound states which host itinerant Bogoliubov 
quasiparticles. These quasiparticles are solid-state real- 
izations of massless Majorana fermions. 

There is currently an intensive search for topologi- 
cal superconductors. In particular, a recently discov- 
ered superconductor Co^E^Ses with T c ~ 3-RT[T3] has 
attracted much attention[TS]. A theoretical study [TT] 
proposed that the strong spin-orbit coupled band struc- 
ture of Cu^E^Sea favors an odd-parity pairing symme- 
try, which leads to a time-reversal-invariant topologi- 
cal superconductor in three dimensions. Subsequently, 
many experimental and theoretical efforts P^-t^O] have 
been made towards understanding superconductivity in 
Cua;Bi2Se3. In a very recent point-contact spectroscopy 
experiment, Sasaki et aL[T] have observed a zero-bias con- 
ductance peak which strongly indicates unconventional 
pairing[2"T]. 

In this Letter, we find a new branch of gapless surface 
Andreev bound states (SABS), in addition to linearly 
dispersing Majorana fermions at k = 0, in the topo- 
logical superconducting phase of Cu2;Bi2Se3 and related 
doped semiconductors. This new branch of SABS is lo- 
cated near the Fermi momentum and is protected by a 
new bulk topological invariant. Moreover, they result 
in unique features in the tunneling spectra which are in 
good agreement with the point-contact spectroscopy ex- 
periment on Cua;Bi2Se3 pQ. We conclude by predicting 
clear signatures of these SABS, which can be tested in 
future tunneling and photoemission experiments at low 
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FIG. 1: a) Side view of a semi-infinite crystal of Bi2Se3. The 
two relevant p z orbitals are shown in the zoom-in view of the 
QL unit cell, b) Bulk and surface bands of the tight-binding 
model for Bi2Se3. [ii and pi 2 denote two chemical potentials 
where the surface states have, respectively, not merged and 
merged into the bulk bands. 



temperatures. 

We start from the k-p Hamiltonian for the band struc- 
ture of Cu x Bi 2 Se3 near r[TT] 



H(k) 



F v zk z Cy + vo~ z (k x Sy — kyS x ). 



(1) 



Here a z = ±1 labels the two Wannier functions which 
are primarily p z orbitals (from Se and Bi atoms) on the 
upper and lower part of the quintuple layer (QL) unit 
cell respectively (see Fig.l). Each orbital has a two-fold 
spin degeneracy labeled by s z = ±1. We note that an 
earlier k-p Hamiltonian[22j violates the mirror symmetry 
of the lattice, and a corrected version[53] is consistent 
with ([I]) . Detailed discussion of the discrepancy is left to 
Supplementary Material [24] . The sign of mv z is a crucial 
quantity which will now be inferred from the existence of 
surface states near k x — k y = in the surface Brillouin 
zone. 

Consider a semi-infinite 0112,612863 crystal occupying 



2 



z < 0, which is naturally cleaved between QLs (see 
Fig.l). The realistic boundary condition corresponding 
to such a termination in the continuum k-p theory is|llj 



a z ip(z = 0) = x/j(z = 0). 



(2) 



This boundary condition reflects the vanishing of the 
electron wavefunction on the bottom layer (er 2 = — 1) 
at z = 0. Solving the differential equation 



Etp = H(k x ,k y ,- 



-id z )ip 



(3) 



subject to (38), we find two branches of mid-gap states 



^±{k x ,k y ,z) = e*/ , (l,0) o .®(l,±ie i *)., (4) 

where I = —v z /m is the decay length, <j) is the azimuthal 
angle of (k x , k y ), and the subscripts a and s denote the 
orbital a z and spin s z basis. For v z m > 0, there are 
no decaying solutions; only when v z m < in Q do 
we obtain surface states decaying in the — z direction. 
The dispersion of these surface states is E±(k x ,k y ) — 



±v^Jk 2 + fc 2 = ±vk, which agree well with the photoe- 

mission data from Cu x Bi2Se3[16 . Thus, the existence of 
surface states on surfaces terminated between QLs estab- 
lishes v z m < in H(k) for Cu x Bi 2 Se 3 [21]. 

Having established that v z m < and v parameterizes 
the linear dispersion of the surface states, we now turn 
to the superconducting state of Cu^t^Ses. Ref.fTTj clas- 
sified four different pairing symmetries compatible with 
short-range pairing interactions, and found that a spin- 
triplet, orbital-singlet, odd-parity pairing symmetry is fa- 
vored when the inter-orbital attraction exceeds the intra- 
orbital one. The mean-field Hamiltonian of this super- 
conducting state is 



Here cj^ — (c^ c^. c k 2^7 ^k 24.) k = c — k 'isy 

four- component electron operators, with the subscript 
1, 2 labeling the two orbitals (Fig. la). In the Bogoliubov- 
de Gennes Hamiltonian H(k), t x and t z are Pauli ma- 
trices in Nambu space, A is the pairing potential, and 
fi > \m\ is the chemical potential in the conduction band. 

The above odd-parity superconducting Cu2;Bi2Se3 is 
fully-gapped in the bulk but has topologically protected 
surface Andreev bound states. To determine the wave- 
function and dispersion of these bound states, we begin 
by solving the BdG Hamiltonian H(k x , k y , —id z ) for the 
SABS at k x — k y = 0. We find a Kramers pair of e = 
eigenstates[2"4"]: 

^fc=o,a(«) = e z - A /^{sin{k F z - 9),sm(k F z)) a 

<8> [(l,-a) a ,isgn(w a )(l,a)J T , a = ±1 (6) 



W(k) 



dk[4,c_k]ft(k) 



Ck 



(H(k) - [L)T Z + A<7 y S z T x . 



(5) 



Here kp = \] \i 2 — m 2 /v z is Fermi momentum in the z 
direction, and 6 is defined by e 10 = (m + i^J fi 2 — m 2 )/ '/z. 
The subscript r denotes a Nambu spinor. The Bo- 
goliubov quasiparticle at k = is defined by 7 a — 
J dz tl>k=o,a(z)(c* (2); c(z))jT =0 . It is straight-forward to 
verify that = j a up to an unimportant overall phase. 
This means that such quasiparticles are two-component 
massless Majorana fermions in 2 + 1 dimensions. 

Having found the SABS wavefunction at e = 0, k = 0, 
we now show that the SABS dispersion crosses e = 
again at finite k, which is one of the main results of this 
paper. We establish this second crossing in two different 
ways: first, by a direct calculation, and second, by a 
topological argument. It will become evident that the 
two approaches yield complementary information. 

In the direct approach, we search for a second crossing 
by asking for which ko > does H(0, fco, —id z )ip = have 
a solution (it suffices to consider k x = 0, k y = ko > 
only, due to rotational invariance). We find that k is 
the nontrivial solution of the algebraic equation [53] 



\x\ 2 + 2sgn{v z )— Re(x) -1 = 0, 
m 



(7) 



where x is defined as 

_ vko — i(A + iEp) 



= , E F = vV-m 2 . (8) 



^/(vk ) 2 + (A + iE F y 

For 0^612863 in the normal state with A = and 
v z m < 0, the above equation has a solution k = ji/v, 
which exactly correspond to the topological insulator sur- 
face states at Fermi energy obtained earlier in Q. With 
superconductivity, topological surface states in the nor- 
mal state turn into SABS, with their location fco and 
wavefunction tj)k ,a perturbed by A: fc ~ £(1 - ^2) 
and ipk ,a acquires particle-hole mixing to first order in 
A. Due to rotational invariance of the fc • p Hamilto- 
nian, the second crossing, hereafter denoted by fco, exists 
along all directions in the xy plane. This leads to a Fermi 
surface of SABS. 

In the topological approach, we first solve for the SABS 
dispersion at small fc and use topological arguments to 
infer its behavior at large fc. Again, we set k x — for 
convenience. Treating the fc y -dependent term in i?Bdc 
as a perturbation, we find the dispersion is linear near 
fc = 0: e a (k) = avk + o(fc 3 ), forming a Majorana cone. 
The velocity v is given by: 



A 2 + sgn(u z )Am 
A 2 + sgn(v z )Am + \i 2 



. . mA 
v ■ sgn{v z )^-. (9) 



In the second equality, we have used the fact A <C |m| < 
fi for weak-coupling superconductors. 

In ([9]), it is important that the SABS velocity v at 
fc = has an opposite sign from the band velocity v 
in the normal state of the doped topological insulator 
Cu x Bi2Se3 (v z m < 0). As we now show, this fact has 
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crucial implications for the SABS dispersion away from 
k = 0: the two branches of SABS t/'fc.i must cross each 
other at e = an odd number of times between T and the 
surface Brillouin zone edge M. The existence of such ad- 
ditional crossings is dictated by a topological invariant we 
call "mirror helicity" , which is a generalization of mirror 
Chcrn number [25] in topological insulators to topologi- 
cal superconductors. To define this invariant, note that 
the crystal structure of Cua;Bi2Se3 has a mirror reflection 
symmetry x — > —x. As a result, the band structure 
is invariant under mirror. However, the pairing poten- 
tial in (|5| changes sign under mirror reflection. So the 
BdG Hamiltonian is invariant under a mirror reflection 
combined with a Zi gauge transformation A — > — A: 



a) 



E./d 



b) 



E/i 



H{k x 



h k 



MH(-k x ,k y ,k z )M- 



(10) 



Here M = Mt z , M = —is x represents mirror reflec- 
tion on electron spin. Because of this generalized mirror 
symmetry, bulk states are grouped into two classes with 
mirror eigenvalues ±z respectively. Each class can have 
a nonzero Chern number n±i. Time reversal symmc- 



The magnitude \n 
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try requires n +i = —r 
determines the number of helical Andreev modes with 
k x = on the edge of yz plane, while the sign de- 
fines a Zi mirror helicity: r\ = sgn(ra + i) = — sgn(n_;). 
The bulk topological invariant rj determines the helic- 
ity of such Andreev modes. For instance, rj < im- 
plies that the mode with mirror eigenvalue —i(+i) moves 
clockwise(anti-clockwise) with respect to +x axis at the 
edge of the yz plane, and its energy-momentum disper- 
sion curve must eventually merge into the E > bulk 
quasiparticle continuum at a large positivc(ncgative) mo- 
mentum. Similar bulk-boundary correspondence applies 
to surface states in topological insulators 25, 26 . 

As we show in Supplementary Material [24], the topo- 
logical superconducting phase of Cu,EBi2Se3 and the un- 
doped topological insulator Bi 2 Se3 have the same mir- 
ror helicity 77, which is determined by the sign of the 
Dirac band velocity v in the bulk. Given the relation 
between 77 and helicity of surface excitations, this implies 
that the SABS in Cu x Bi2Se3 must have the same helic- 
ity as surface states in Bi 2 Se3. On the other hand, the 
SABS velocity v at k = has an opposite sign from the 
Dirac band v. To reconcile this fact with the helicity re- 
quirement, the two SABS branches ipk,a — which are mir- 
ror eigenstates with eigenvalues M = ia — must become 
twisted and switch places before merging into the bulk. 
This necessarily results in an odd number of crossings 
between f and M. 

The above topological argument reveals the robustness 
of gapless SABS at the second crossing in the k-p regime 
and beyond. In the k ■ p regime, the surface states at k 
and — k have opposite mirror eigenvalues (or spins) due 
to their helical nature, whereas the pairing symmetry 
A only pairs states with the same mirror eigenvalues. 
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FIG. 2: SABS dispersion for the tight-binding model in which 
a) m = —0.3 < 0, fxi = 0.6 and b) m — —0.3 < 0, [i% = 
1; The mirror eigenvalues are displayed near each branch of 
SABS. The SABS twist with a second crossing near Fermi 
momentum, as also observed in Ref. 20 . The arrow denotes 
where the dispersion has zero slope, resulting in a Van Hove 
singularity in the density of states. 



This symmetry incompatibility makes the surface states 
remain gapless in the topological superconducting phase 
|27) . Moreover, the topological argument demonstrates 
that the second crossing is topologically protected by the 
mirror helicity invariant in the bulk, as long as v/v < at 
k = 0. As a result, the second crossing remains in a much 
larger energy range, even when higher order corrections 
to the k ■ p Hamiltonian become important, as shown 
below. In particular, we emphasize that the existence 
of the second crossing is independent of whether surface 
states are separated from the bulk at the Fermi energy. 

To gain more insight into these twisted SABS and to 
calculate their local density of states, we explicitly obtain 
its dispersion in the entire surface Brillouin zone. For this 
purpose, we construct a two-orbital tight-binding model 
in the rhombohedral lattice shown in Fig.l and calculate 
the SABS dispersion numerically. Details of our tight- 
binding model and its distinction from previous models [1 , 
|2"0] are described in the Supplementary Material [24). 

Here we would like to note the following aspects of 
our model. The normal state tight-binding model is con- 
structed to reproduce both the k ■ p Hamiltonian (1) of 
011^612863 in the small k limit and the boundary condi- 



tion ( 38 ) in the continuum theory. The bulk and surface 
bands of the normal state tight-binding model are dis- 
played in Figure lb; at chemical potential (Xi, the Fermi 
momentum is relatively small and terms higher order 
than k are negligible, whereas at /J2, these higher order 
terms cause deviation from the k ■ p Hamiltonian. 

Upon adding odd-parity superconductivity pairing to 
the model, we obtain the SABS dispersion (Fig. 2). A 
branch of linearly dispersing Majorana fermions is found 
at k = 0, which signifies a three-dimensional topologi- 
cal superconductor. In addition, the bands of Andreev 
bound states in the surface Brillouin zone are twisted: 
they connect the Majorana fermion at k = with the 
second crossing near Fermi momentum. Such behavior 
was independently found by Hao and Lee [5111 121], ancl its 
topological origin is revealed by our analytical calcula- 
tions and arguments. 
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For a given branch (M = ±i) of SABS, its particle-hole 
character evolves as a function of momentum from hav- 
ing an equal amount of particle and hole (charge neutral) 
at k = to being exclusively hole or particle (charged) at 
large k. At chemical potential /ii, the SABS near the sec- 
ond crossing can be identified with nearly unpaired sur- 
face states in the normal state, which show up twice — as 
particle and hole — in the BdG spectrum. However, even 
when these surface states have merged into the bulk, the 
SABS still has the second crossing, as required by the 
mirror helicity. This is shown in Fig. 2b, at chemical 
potential The resulting gapless SABS near the sec- 
ond crossing has substantially more particle-hole mixing 
than the first case and is unrelated to surface states in the 
normal state. Such SABS defy a quasi-classical descrip- 
tion and represent a new type of Andreev bound states 
which arises from the interplay between nontrivial band 
structure and unconventional superconductivity. 

Finally, we relate our findings of SABS in Cuj;Bi2Se3 
to the recent point-contact spectroscopy experiment pQ , 
in which a zero-bias differential conductance peak along 
with a dip near the superconducting gap edge was ob- 
served below 1.2K and attributed to SABS. To com- 
pare with this experiment, we calculate the local tun- 
neling density of states (LDOS) as a function of energy 
for m//i 2 = 0.3 — roughly the value found in ARPESfTB], 
The resulting LDOS at zero and finite temperatures are 
shown in Fig. 3. The finite temperature LDOS from 
T = 0.05A to T = 0.2A agrees with the experimen- 
tally observed differential conductance peaks as well as 
the dips with the slight asymmetry between positive and 
negative voltages. Both features along with the absence 
of coherence peaks contrast sharply with the tunneling 
spectrum of an s-wave superconductor. 

In addition to comparison with the experiment, we 
make the following predictions stemming from the zero 
temperature LDOS in Figure 3a. Here the two peaks 
arise from Van Hove singularities at the particular en- 
ergy near E = where the SABS bands have zero slope, 
indicated by the arrow in Fig. 2b. Furthermore, the 
significant asymmetry in the height of these two peaks 
reflects the fact that the SABS at the turning point is 
primarily of hole type, as noted earlier. The energy of 
these two peaks and the magnitude of their asymmetry 
depends somewhat on details of band structure. How- 
ever, the existence of two peaks only depends on there 
being a turning point in the SABS dispersion, which is 
guaranteed by the existence of a second crossing in a wide 
regime of chemical potentials. Hence, we predict that 
for relatively clean surfaces the zero-bias conductance 
peak in the tunneling spectra will split into two asym- 
metric peaks at even lower temperatures. Such peaks 
will be an unambiguous signature of Majorana fermions 
smoothly turning into normal surface electrons. Further- 
more, the SABS dispersion we predict in Fig. 2 can be 
directly tested in future ARPES experiments. 




FIG. 3: Tunneling local density of states (arbitrary units) 
at a) T — and b) finite temperature. In both cases, the 
chemical potential is fi% = 1. 



While the main focus of this Letter is Cu x Bi2Se3, 
we end by discussing the implications of our findings 
for superconducting doped semiconductors with similar 
band structures. Candidates include Bi 2 Te3[3Tj under- 
pressure, TlBiTe 2 [32], PbTe[33], SnTe[33, and GeTe[35]. 
Provided that the material is inversion symmetric and 
its Fermi surface is centered at time-reversal-invariant 
momenta, the Dirac-type relativistic k ■ p Hamiltonian 
([!]) describes their band structures [28]. Moreover, if the 
pairing symmetry is odd under spatial inversion and fully 
gapped, the system is (almost) guaranteed to be a topo- 
logical superconductor according to our criterion |1U 130] . 
Our work is also relevant to noncentrosymmetric super- 
conductors such as YPtBi 36J, if their pairing symmetries 
have dominant odd-parity components. 

As a final point which captures the essence of this work, 
we compare and contrast SABS in doped superconduct- 
ing topological insulators with normal insulators, which 
differ by a band inversion (v z m < versus v z m > 0). 
In both, the Majorana fermion SABS exist at k = as 
shown in ([6j [9|) . However, the SABS in doped normal in- 
sulators do not necessarily have the second crossing near 
Fermi momentum[24 . This can be understood from our 
mirror helicity argument, with the difference being that 
v/v > for n 2 m > (see Eq.Q). In this sense, the 
new type of surface Andreev bound state and its phe- 
nomenological consequences are the unique offspring of 
both nontrivial band structure and odd-parity topologi- 
cal superconductivity. 

Note: Two recent studies [TJ [2U] calculated the surface 
spectral function numerically in Cu x Bi2Se3 tight-binding 
models. The second crossing of SABS was independently 
found in Ref. |20j . We also learned of another point- 
contact measurement on Cu^Bi9Se^ |37j. 
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SUPPLEMENTARY MATERIAL 
I. SABS Wavefunction and Second Crossing 

First, we derive in detail the wavefunction (6) in the 
main text from the BdG Hamiltonian 'H(k). A Kramers 
pair of zero-energy eigenstates tpk=a,a=± (z) with mirror 
eigenvalues M — i ■ a is expected from the topology and 
symmetry of %(k|| = 0). For a given mirror eigenstate, 
s x is locked to r z by the identity s x t z = —a, so that 
ipk=o,a(z) satisfies a reduced 4-component equation: 



[(ma x 



Jz a y d z - y)r z + Aa y T x }iJj k=0ya (z) = 0. (11) 



This can be further simplified by multiplying both sides 
by t z : 

[ma x - iv z a y d z - fi + iAa y Ty]tJj k=0 ^ a (z) = 0. 

It is evident that ipk=o,a is an eigenstate of r y . The 
corresponding eigenvalue is given by sgn(v z ) in order to 
have a decaying solution. Eq.( 11 ) then reduces to a two- 



component equation in orbital space, which has two in- 
dependent solutions: 

£±(*0 = (h e ±l6 )a ■ e (±^+ A /M) z . (12) 



9 is defined by e 10 = (m + i\J 'fi 2 — m 2 )/fi. Choosing a 
suitable linear combination of £ + and to satisfy the 
boundary condition (2) in the main text, we obtain the 
wavefunction of SABS, which is reproduced here for the 
reader's convenience: 



lpk=0,a(z) 



(sin(/c^z — 9), sm(kpz)) a 



® [(1, -a) s ,isgn{v z )(l,a) 



±1. 



Next we solve for the location of the SABS second 
crossing. For convenience, we look for a zero-energy so- 
lution ip{z) at k x = 0, k y = ko with mirror eigenvalue +i 
(i.e., s x t z — —1). ip satisfies 

[{m<j x - iv z <T y d z - p)t z + vk a z + Aa y T x ] ip(z) = 0(13) 

Recall that v z m < for a doped topological insulator. 
Without loss of generality, here we choose m < 0, v z > 0. 



By multiplying Eq.( 13 1 by ia y T z , the zero-energy solution 
satisfies 

[ma z + v z d z - iiiOy - vk a x r z - Ar y ] ip(z) = 0. (14) 
We write the wavefunction ip(z) as 

$(z) = e lXa ^ 2 (j){z), A € C (15) 
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where cosA = ulEp and sinA = —im/Ep, Ep 
y/ H 2 — m 2 . Eq.(13l then becomes 



[v z d z - iEpOy - vk a x r z - Ar y ] , 



(16) 



Note that Eq.(16) now commutes with cr y T y , which be- 



comes a constant labeled by r. The reduced equation for 

(j) T {z) is 

[v z d z - (rA + iEp)<7 y - vk <J x ] cj>{z) = 0. (17) 

The solution takes the form <j>{z) = e Kz £. First consider 
t = 1. From Eq. (17), we have 



K± = W^±W±viZ s±E , /v ,. (18) 

Corresponding eigenvectors are given by 

C± = (as±, 1), x± = (vk - i(A + iE F ))/{v z K ± ) (19) 

To get a decaying solution, we must have Re(iT) > 0. 
Hence, we must choose K + and thus £+. We now rewrite 
the complete wavefunction with both orbital and Nambu 
components (spin is locked by s x t z — — 1 and not shown 
explicitly) : 

y- X + _ • ■ ,1 X + + i f-, ■ ■ „ 

4+ = — 2 — (M,«,-1)H o — (1> -*,-*>-!) 

= (a; + , 1,1,-1+) (20) 

Note that the equation for r = —1 is equivalent to the 
complex conjugate of that for r = 1. Therefore if we 
choose K+ = K and x + = x for t = 1, we must choose 
K* and x* for r = — 1. The corresponding wavefunction 



e = (x*,i,i,- x *) 



(21) 



It follows from Eqn. (15) that 



ip(r = 1) = e lXa */ 2 (x,lA,-x) 
V(r = -1) = e lA ^/ 2 (x*,l,-l,x*) (22) 

Up to normalization, the most general form of ip(z) sat- 
isfying the boundary condition ( 38 ) is 



V>(0) = (i,o, AO), 



(23) 



where ^4 is some constant. Hence, for a nontrivial solu- 
tion to exist, the determinant of the 2x2 matrix made 
from the second and fourth component of ip(r = 1) and 
ip(r = —1) must be zero. This condition is simplified to 
an algebraic equation 

= 2Rc(.t) + ^(-1 + |x| 2 ) (24) 
Ep 

which is the result cited in the main text. Our previ- 
ous solution at k — ^ corresponds to x = ±i, which 



satisfies the above condition. Another simple limit is 
the normal state with A = 0. In this case, the second 
crossing is simply located at the momentum where the 
topological insuator surface states cross the chemical po- 
tential, namely, fco = fJ-/v. We can check that for this 
case, x (fi + Ep)/{—m) indeed satisfies Eq. ( p4[ ). Now 
we solve for fco to first order in A. Temporarily absorb- 
ing v into fco and expanding x to second order in A, we 
obtain 



Re(x) = 



E h 



fc 



(-2E F k Q - fc 2 )A 2 



Im(x) 



2(E F - k Q ) 2 (E F + k )^/-E 2 F + k 2 
k A 



m 
~E~p 



{E F - fc ) v^F+^o 
-2m 2mfc A 2 



fc [E F -k ) 3 (E F + k ) 



From Eq. ( 24 ) , we then extract the leading order to cor- 



rection to fcn 



fco 



(1 



2m 2 



The corresponding x at fco is given 
Re(a:) = 
Im(a;) = — 



m 

-fx + E F 



A 2 n(n + E F ) 
2m 3 (-fi + E F ) 



A/i 



m(— /U + E F ) 



(25) 

(26) 
(27) 



We conclude this section by calculating the ratio of the 
particle (r = 1) and hole (r = —1) components of the 
s x t z — — 1 wavefunction ip(z) at the second crossing and 
at z = 0. This wavefunction is some linear combination 
c\ij)(T = 1) + c<i{t = —1) with vanishing second and 
fourth components (to satisfy the boundary condition). 
Hence, we find 



c 2 — (cos(A/2) + ix sin A/2) 
ci cos A/2 + ix* sin A/2 

The hole/particle ratio is 

cos A/2(ci — c 2 ) — i sin \/2(c\x — c 2 x*) 

T — 

cos(A/2)(cix + c 2 x*) + isin(A/2)(ci + c 2 ) 



(28) 



(29) 



Using the fact that the second and fourth components 
vanish, which is equivalent to 

cos A/2(ci + c 2 ) + zsin \/2(cix + c 2 x*) = (30) 
— cos \j2(c\x — c 2 x*) + zsin A/2(ci — c 2 ) = 0, (31) 



we get 



ci - c 2 



(ci +c 2 )(icot(A/2)) 



1 + i(tanA/2)Re(ir) 
ilm(x) 



(32) 
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FIG. 4: Phase diagram of fully-gapped odd-parity supercon- 
ductivity in doped semiconductors as a function of band gap 
m and pairing potential A, showing three gapped phases: 
band insulator, topological insulator and topological super- 
conductor. They are topologically distinguished by the mirror 
Chcrn numbers n+i. 



Recalling that cos A = fi/Ep and sin A = —im/Ep, we 
have 



tan(A/2) 



(i + m — Ep 
i(e iX + 1) ~ fi + m + E F 



(33) 



The hole/particle ratio at the second crossing is thus 
A(n + E F )(fi + m - E F ) 



r = 



(34) 



2im 2 (fi + m + E F ) 
which is first order in A. 

II. Mirror Helicity 

Here we show that the topological insulator and topo- 
logical superconductor phases have the same mirror he- 
licity. We deduce this fact from the phase transition 
between topological insulators and topological supercon- 
ductors. 

The BdG Hamiltonian (5) in the main text exhibits 
three topologically distinct gapped phases as a function 
of the band gap, pairing potential and doping. At zero 
doping (/i = 0) and in the absence of superconductivity 
(A = 0), the system is either an normal band insulator or 
a topological insulator, depending on the sign of m. At fi- 
nite electron doping, the chemical potential lies inside the 
conduction band: fi > 0. When the odd-parity pairing 
A occurs in such a doped normal insulator or topological 
insulator, the system becomes a fully gapped topological 
superconductor. For the sake of our argument, we note 
that the topological superconductor phase is adiabati- 
cally connected to the \i = and A > \m\ limit. Fig. 4 
shows the three phases in the fi = phase diagram as 
a function of m and A. The phase transition between 
topological superconductors and normal/topological in- 
sulators occurs at A = ±ra. 

Recall from the main text that due to mirror symme- 
try, each phase has a mirror Chern number n+j = n_j 



displayed in Fig. 4. Using n + i = for the normal in- 
sulator as a reference, we can obtain the mirror Chcrn 
number for the topological insulator and topological su- 
perconductor by calculating the change of n +i across the 
phase transition to the normal insulator. Due to the 
double counting of particles and holes, the mirror Chern 
number of a band insulator defined in Nambu space is 
always an even integer twice the value of that defined 
previously for insulators in Ref.[25]. As a result, a direct 
transition from topological insulator to band insulator at 
A = changes n +i by two. For A^O, this transition is 
split into two transitions with an intermediate topologi- 
cal superconductor phase, so that each transition changes 
n + i by one. Therefore we have 



n +i (IT) = 2n+,(TSC). 



(35) 



Recall that mirror helicity is defined as 77 = sgn(n+j). 
Hence, the topological insulator and topological super- 
conductor phase have the same mirror helicity. 



III. Tight-binding Model 

Here we present the details of our tight-binding model. 
This model is defined on the rhombohcdral lattice with 
a bilayer unit cell shown in Fig.l. The Hamiltonian H = 
H a +Hi2 + H soc +H' 12 consists of the following four terms. 



<ij> 



describes nearest neighbor hopping within the same layer. 



H 12 = tlC la c j a 

<i£l,j£2> 



*2ct a c j / Q (36) 

<i£l,j'£2> 



describes hopping between two adjacent layers within a 
QL (ti) and on two neighboring QLs (£2). to, t\ and t 2 
are spin-independent. In addition, the two orbitals in the 
upper and lower part of the unit cell (Fig. la) experience 
local electric fields along the ±z direction, which give rise 
to the following Rashba-type spin-orbital associated with 
intra-layer hopping: 



( E 

<ij>ei 



E ) 

<ij>E2 



iX + _ 



pCjp ■ (z x a tJ ), 



where = ^Cijk(Rj — Rfe) denote the vectors 

joining nearest neighbors within a layer, and R.1,2,3 
are the Bravais lattice vectors. The last term H' 12 
(which plays a minor role) describes inter-layer sec- 
ond nearest neighbor (i 3 ) hopping within a QL: H' 12 = 

^2<<i£l,j£2>> ^3 c ia c ja + h.C. 

We emphasize that our tight-binding Hamiltonian H, 
by construction, satisfies the symmetries of the Bi 2 Se 3 
crystal structure. Its point group D%d has three inde- 
pendent symmetry operations: inversion P, three-fold 



rotation C3 around the z axis, and reflection M about 
the x axis. These operations act on the orbital and spin 
degrees of freedom as follows: P interchanges the two 
orbitals (see Fig. la), C3 rotates the electron spin s x and 
s y , and M flips s z and s y , but not s x (Recall that spin 
is a pseudovector) . Therefore, these operations are rep- 
resented by P — <r x ,Cz — exp(—i^^-s z ), M = —is x . 

The above tight-binding model captures the essential 
features of the band structure Bi2Se3 near the T point. 
(We caution the reader that our tight-binding model does 
not aim to describe the band structure of Cu-rE^Sea in 
the entire Brillouin zone. Such a task requires a real- 
istic band structure modeling beyond the scope of this 
work.) First, the Bloch Hamiltonian H(k) reduces to 
the k ■ p Hamiltonian (1) in the main text as follows: 
m = 3(ti + t 2 + t 3 ), v z = 3^2C, and v — |Aa 2 , where 
a = jay I and c = || (Ri + R2 + R3) | . Second, our model 
is able to reproduce the Dirac surface states (Fig. lb). To 
understand this, we note that at k x = k y = 0, the spin- 
orbit term H soc vanishes. The resulting one-dimensional 
system corresponding to H(k x = k v = 0) is equivalent 
to the well-known Su-Heeger-Schrieffer model for poly- 
acetylene, which has a similar two-site unit cell. In both 
systems, the hopping between neighboring sites within a 
unit cell is different from that between two unit cells. As 
a result, when such a one-dimensional system is termi- 
nated on a "strong bond", zero-dimensional end states 
appear within the band gap and are spin degenerate. 
In contrast, when the system is terminated on a "weak 
bond", end states are absent. In the context of Bi2Se3, 
strong bond correspond to termination between two QLs, 
and weak bond correspond to termination within a QL. 
In the former case, the end states at k. x = k y = dis- 
perse and become spin-split as a function of k x and ky, 
due to the fc-linear spin-orbit term H soc . As a result, 
they constitute the two-dimensional Dirac surface band 
of Bi 2 Se 3 . In the latter case, the end states are absent 
at k x = k y = 0. Instead, the surface state Dirac points 
of Bi 2 Se3 are located at three M points [251 I3"8"] (which 
cannot and should not be accessed by k ■ p Hamiltonian 
near T). It will be interesting to experimentally verify 
such a drastic dependence of surface band structure on 
surface terminations. 

To capture the effect of two different surface termina- 
tions within a continuum theory, we choose the boundary 
condition correspondingly. The boundary condition for 
termination between two QLs (strong bond) is 

a z ip(z = 0) = iP(z = 0). (37) 

This reflects the vanishing of the o~ z = — 1 component 
of the wavefunction at z = (the outmost site corre- 
sponds to a z = 1). Instead, the boundary condition for 
termination with a QL (weak bond) is 

<t z tP(z = 0) = -tP(z = 0). (38) 




FIG. 5: SABS dispersion for the tight-binding model in which 
a) m = 0.3 > 0, /ii = 0.6 corresponds to a doped BI; b) m = 
0, jUi = 0.6 corresponds to a doped zero-gap semiconductor. 

As we have shown in the main text, for v z m < Dirac 
surface states exist in the continuum theory for the first 
termination, but not for the second. This correctly re- 
produces the experimental phenomenology. 

To include superconductivity, we add the following 
odd-parity pairing term in the Hamiltonian: 

h M f = h+ j + c lsV) + h - c - 

<i£l,j£2> 

The parameters we used are A = 0.03, to = — 0.1, £1 = 
-1,* 2 = 0.5, t 3 = 0.6, a = l,c = 1,A = 0.5, and fi = 0.6 
(above the normal state surface Dirac point) for Figure 
2a and [i — 1 (above the Dirac point) for Figure 2b. The 
slab size was 320 unit cells. We note that v z oc i 2 > 
actually corresponds to v z < in the k ■ p Hamiltonian 
above because our simulated crystal is oriented in the 
opposite z direction relative to the k ■ p definition. 

For completeness, we calculate the SABS dispersion for 
a doped band insulator (mv z > 0), in which the second 
crossing does not exist because v/v > (Fig. 5a). The 
dispersion for the critical case (to = 0) is displayed in 
Fig. 5b. 

IV. Relation to Previous Works 

In a recent work, Hao and Lee [50] calculated the 
surface spectral function in tight-binding models for 
011^612803 with the four possible pairing symmetries!!!), 
including the fully-gapped odd-parity pairing studied in 
this work. They used two tight-binding models which 
are lattice regularizations of two k ■ p Hamiltonians (I 
and II). However, both these Hamiltonians violate the 
mirror symmetry. Model II is quoted from the incorrect 
k ■ p Hamiltonian of Ref. 22J: their terms k z a x s z as well 
as <J x (k x s x + k y Sy), in the basis they specify, violates the 
mirror symmetry M. A corrected version[23 is identical 
to our k ■ p Hamiltonian (1) after interchanging o~ x and 
a z (corresponding to a change of basis for the orbitals) . 
Model I is claimed to be quoted from Ref.fTT] (the one 
we use here). However, the term o~ z (k x s y — k y s x ) is mis- 
takenly replaced by a z (k x s x + k y s y ). 

Nonetheless, if one forgoes the definition of s x y as op- 
erators corresponding to spin along the x and y direc- 
tions in real space, then their Model I corresponds to our 
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k-p Hamiltonian after a unitary spin rotation exp(— ijs z ) 
(without affecting the odd-parity pairing term Aa y s z T x ). 
Hence, they aiso found that Majorana fermion Andreev 
bound states at k — connect to the Dirac surface 
states near Fermi momentum. They attributed the sec- 
ond crossing to the fact that Dirac surface states remain 
gapless in the odd-parity superconducting state, and con- 
cluded that it disappears if the surface states merge into 
the bulk. In contrast, our work revealed the topological 
origin of the twisted surface Andreev bound states: as 
long as v/v < 0, they are protected by mirror symmetry 
and exist independent of whether Dirac surface states 
appear at Fermi energy. 



where > is the energy cost of creating a quasi- 
particle excitation, u and v are the particle and hole 
components of the positive-energy eigenstates of BdG 



Hamiltonian, respectively. To derive (42), one must keep 
in mind that adding (removing) a quasi-particle always 
increases(decreases) the energy of the system. Because 
the hole component of a E > eigenstate is related to 
the particle component of its partner at — E by the inher- 



ent particle- hole symmetry in BdG formalism, ( 42 1 can 
be simplified to 

A~(uj, k) = \u k \ 2 n F (uj)S(u! - £ fc ), 

A+{oj,k) = \u k \ 2 (l-n F (oj))6(u-£ k ), (42) 



V. Finite Temperature Differential Conductance 

Finally, we elaborate on how we attained the differen- 
tial conductance plots in the main text. Consider two 
systems separated by an insulating barrier. Then the 
tunneling current is proportional to the transition rate 
given by Fermi's golden rule: 

Ioc J deAf(e + eV)A^(e) - A± {e + eV) A% (e) (39) 

where A + (e) is the probability of adding a particle and 
changing the system's energy by e (positive or negative), 
and A~ (e) is the probability of removing a particle and 
changing the system's energy by — e. 1 and 2 denote the 
two sides of the barrier. 

For free electron systems, A ± is given by the density 
of states weighted by the Fermi-Dirac distribution 



A±(e) 



dkA ± {e 1 k) 



A-(e,k) = n F (e)5(e-£ k ) 
A+(e,k) = (l-n F (e))S(e-£ k ) 



(40) 



where n F (e) is the Fermi-Dirac distribution function 
l/(e e / T + 1). For convenience, hereafter both e and ^ k 
are measured with respect to chemical potential. 
For a BCS superconductor, A^ is modified: 

A-(oj,k) = \u k \ 2 n F (co)S(uj - l&l), lu > 

M 2 (i-n F (MMM-|ai), ^<o 

A+(u,k) = \u k \ 2 (l-n F (u>))6(u]-\$ k \), lj>0 

|» fc |V(M)«(M-|e*|), u;<0 (41) 



where we have used 1 — n F (— uj) = n F {uj). Here ui and 
£ k can be both positive and negative. Written in this 
form, A for a superconductor is similar to a normal 
metal, except it has prefactor u k . When superconduc- 
tivity vanishes, u k — l,v k — for k > k F ,cj > and 
k < k F ,oj < 0, whereas u k = 0,v k = 1 for k < k F ,cu > 



and k > k F ,uj < 0. In this limit, (42) reduces to the free 



fermion case ( 40 ) 



In our simulation, \u k \ 2 and \v k \ 2 were obtained from 
the r = 1 and r = — 1 components of the surface Green's 
function, summed over spin and for the a z = +1 orbital 
at z = only, in accordance with our boundary condi- 
tion. The surface Green's function was computed using 
a recursive algorithm [39] , allowing us to use a very large 
slab size (2 layers). 

Substituting A into the expression for tunneling cur- 
rent and assuming that the density of states of the normal 
metal is constant, we obtain 



dep N p s {w)[n F (uj - eV) - n F (u)}, (43) 



Differentiating / with respect to V gives the differential 
conductance 



dl/dV oc / deps{w){dn F / dui)\ u - e v 

J — oo 

p s (u) = [ dk\u k \ 2 S(uj-^ k ). (44) 



